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Abstract 

Absence epilepsy is believed to be associated with the abnormal interactions between the cerebral cortex and thalannus. 
Besides the direct coupling, anatomical evidence indicates that the cerebral cortex and thalamus also communicate 
indirectly through an important intermediate bridge-basal ganglia. It has been thus postulated that the basal ganglia might 
play key roles in the modulation of absence seizures, but the relevant biophysical mechanisms are still not completely 
established. Using a biophysically based model, we demonstrate here that the typical absence seizure activities can be 
controlled and modulated by the direct GABAergic projections from the substantia nigra pars reticulata (SNr) to either the 
thalamic reticular nucleus (TRN) or the specific relay nuclei (SRN) of thalamus, through different biophysical mechanisms. 
Under certain conditions, these two types of seizure control are observed to coexist in the same network. More importantly, 
due to the competition between the inhibitory SNr-TRN and SNr-SRN pathways, we find that both decreasing and 
increasing the activation of SNr neurons from the normal level may considerably suppress the generation of spike-and-slow 
wave discharges in the coexistence region. Overall, these results highlight the bidirectional functional roles of basal ganglia 
in controlling and modulating absence seizures, and might provide novel insights into the therapeutic treatments of this 
brain disorder. 
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Introduction 

Absence epilepsy is a generalized non-convulsive seizure disorder 
of the brain, mainly occurring in the childhood years [1]. A typical 
attack of absence seizures is characterized by a brief loss of 
consciousness that starts and terminates abrupdy, and meanwhile 
an electrophysiological hallmark, i.e. the bilaterally synchronous 
spike and wave discharges (SWDs) with a slow frequency at 
approximately 2-4 Hz, can be observed on the electroencephalo- 
gram (EEG) of patients [1,2]. There is a broad consensus that the 
generation of SWDs during absence seizures is due to the abnormal 
interactions between cerebral cortex and thalamus, which together 
form the so-called corticothalamic system. The direct evidence in 
support of this view is based on simultaneous recordings of cortex 
and thalamus from both rodent animal models and clinical patients 
[3-5] . Recent computational modelling studies on this prominent 
brain disorder also approved the above viewpoint and provided 
more deep insights into the possible generation mechanism of 
SWDs in the corticothalamic system [6-13]. 

The basal ganglia comprise a group of interconnected 
subcortical nucleus and, as a whole, represent one fundamental 
processing unit of the brain. It has been reported that the basal 



ganglia are highly associated with a variety of brain functions and 
diseases, such as cognitive [14], emotional functions [15], motor 
control [16], Parkinson's disease [17,18], and epilepsy [19,20]. 
Anatomically, the basal ganglia receive multiple projections from 
both the cerebral cortex and thalamus, and in turn send both 
direct and indirect output projections to the thalamus. These 
connections enable the activities of the basal ganglia to influence 
the dynamics of the corticothalamic system. Therefore, it is 
naturally expected that the basal ganglia may provide an active 
role in mediating between seizure and non-seizure states for 
absence epileptic patients. Such hypothesis has been confirmed by 
both previous animal experiments [19,21-23] and recent human 
neuroimage data [20,24,25]. Nevertheless, due to the complicated 
interactions between basal ganglia and thalamus, the underlying 
neural mechanisms on how the basal ganglia control the absence 
seizure activities are still remain unclear. 

From the anatomical perspective, the substantia nigra pars 
reticulata (SNr) is one of the major output nucleus of the basal 
ganglia to thalamus. Previous experimental studies using various 
rodent animal models have demonstrated that suitable changes 
in the firing of SNr neurons can modulate the occurrence of 
absence seizures [21-23,26]. Specifically, it has been found that 
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Author Summary 

Epilepsy is a general term for conditions with recurring 
seizures. Absence seizures are one of several kinds of 
seizures, which are characterized by typical 2-4 Hz spike- 
and-slow wave discharges (SWDs). There is accumulating 
evidence that absence seizures are due to abnormal 
interactions between cerebral cortex and thalamus, and 
the basal ganglia may take part in controlling such brain 
disease via the indirect basal ganglia-thalamic pathway 
relaying at superior colliculus. Actually, the basal ganglia 
not only send indirect signals to thalamus, but also 
communicate with several key nuclei of thalamus through 
multiple direct GABAergic projections. Nevertheless, 
whether and how these direct pathways regulate absence 
seizure activities are still remain unknown. By computa- 
tional modelling, we predicted that two direct inhibitory 
basal ganglia-thalamic pathways emitting from the sub- 
stantia nigra pars reticulata may also participate in the 
control of absence seizures. Furthermore, we showed that 
these two types of seizure control can coexist in the same 
network, and depending on the instant network state, 
both lowing and increasing the activation of SNr neurons 
may inhibit the SWDs due to the existence of competition. 
Our findings emphasize the bidirectional modulation 
effects of basal ganglia on absence seizures, and might 
have physiological implications on the treatment of 
absence epilepsy. 



pharmacological inactivation of the SNr by injecting y-aniinobu- 
tyric acids (GABA) agonists or glutamate antagonists suppresses 
absence seizures [21,22]. Such antiepUeptic effect was supposed to 
be attributed to the overall inhibitory effect of the indirect pathway 
from the SNr to thalamic reticular nucleus (TRN) relaying at 
superior colliculus [21,22]. In addition to this indirect inhibitory 
pathway, it is known that the SNr also contains GABAergic 
neurons directly projecting to the TRN and specific relay 
nuclei (SRN) of thalamus [27,28]. Theoretically, changing the 
activation level of SNr may also significantly impact the firing 
activities of SRN and TRN neurons [28,29]. This contribution 
might further interrupt the occurrence of SWDs in the cortico- 
tlialamic system, thus providing an alternative mechanism to 
regulate typical absence seizure activities. To our knowledge, 
however, so far the precise roles of these direct basal ganglia- 
thalamic pathways in controlling absence seizures are not 
completely established. 

To address this question, we develop a reahstic mean-field 
model for the basal ganglia-corticothalamic (BGCT) network in 
the present study. Using various dynamic analysis techniques, we 
show that the absence seizures are controlled and modulated 
either by the isolated SNr-TRN pathway or the isolated SNr-SRN 
pathway. Under suitable conditions, these two types of modula- 
tions are observed to coexist in the same network. Importantly, 
in this coexist region, both low and high activation levels of 
SNr neurons can suppress the occurrence of SWDs due to the 
competition between these two direct inhibitory basal ganglia- 
thalamic pathways. These findings clearly outline a bidirectional 
control of absence seizures by the basal ganglia, which is a 
novel phenomenon that has never been identified both in 
previous experimental and modelling studies. Our results, on the 
one hand, further improve the understanding of the significant role 
of basal ganglia in controlling absence seizure activities, and on the 
other hand, provide testable hypotheses for future experimental 
studies. 



Materials and Methods 

Model 

We buUd a biophysically based model that describes the 
population dynamics of the BGCT network to investigate the 
possible roles of basal ganglia in the control of absence seizures. 
The network framework of this model is inspired by recent 
modelling studies on Parkinson's disease [30,31], which is 
shown schematically in Fig. 1. The network totally includes nine 
neural populations, which are indicated as follows: e = excitatory 
pyramidal neurons (EPN); / = inhibitory interneurons (UN); 
)- = TRN; 5 = SRN; (i|= striatal Dl neurons; (i2 = striatal D2 
neurons; ^i=SNr; /i2= globus paUidus external (GPe) segment; 
C = subthalamic nucleus (STN). Similar to other modelling studies 
[30-32], we do not model the globus pallidus internal (GPi) 
segment independently but consider SNr and GPi as a single 
structure in the present study, because they are reported to have 
closely related inputs and outputs, as well as similarities in cytology 
and function. Three types of neural projections are contained in 
the BGCT network. For sake of clarity, we employ different line 
types and heads to distinguish them (see Fig. 1). The red lines with 
arrow heads denote the excitatory projections mediated by 
glutamate, whereas the blue sohd and dashed lines with round 
heads represent the inhibitory projections mediated by GABAa 
and GABAb, respectively. It should be noted that in the present 
study the connections among different neural populations are 
mainly inspired by previous modelling studies [30,31]. Addition- 
ally, we also add the connection sending from SNr to TRN in our 
model, because recent anatomical findings have provided evidence 
that the SNr also contains GABAergic neurons directiy projecting 
to the TRN [27-29]. 

The dynamics of neural populations are characterized by the 
mean-field model [9,11,33-35], which was proposed to study 
the macroscopic dynamics of neural populations in a simple 
yet efficient way. The first component of the mean-field model 
describes the average response of populations of neurons to 
changes in cell body potential. For each neural population, the 
relationship between the mean firing rate Qu and its correspond- 
ing mean membrane potential Va satisfies an increasing sigmoid 
function, given by 



Qa{r,t) = F[Va(r,t) 



1 +exp 



n {Va{r,t)-6„) 



n/3 



where aeA = {e,i,r,s,d\,d2,p\,P2X} indicate different neural pop- 
ulations, 2™"^ denotes the maximum firing rate, r represents the 
spatial position, is the mean firing threshold, and a is the 
threshold variability of firing rate. If Va exceeds the threshold 0a, 
the neural population fnes action potentials with an average firing 
rate Qa. It should be noted that the sigmoid shape of Qa is 
physiologically crucial for this model, ensuring that the average 
firing rate cannot exceed the maximum firing rate Q^"^ . The 
changes of the average membrane potential Va at the position r, 
under incoming postsynaptic potentials from other neurons, are 
modeled as [9,33-36] 



D,,iVa(t,t)-. 



(2) 



(3) 
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Glutamate 

GABA^ 
GABA, 



Figure 1. Framework of the basal ganglia-corticothalamic network. Neural populations include e = excitatory pyramidal neurons (EPN); 
= inhibitory interneurons (UN); )■ = thalamic reticular nucleus (TRN); 4 = specific relay nuclei (SRN); d\ = striatal D1 neurons; (/2 = striatal D2 neurons; 
p\ = substantia nigra pars reticulata (SNr); = globus pallidus external (GPe) segment; C = subthalamic nucleus (STN). Note that we do not model the 
globus pallidus internal (GPi) segment independently but consider SNr and GPi as a single structure in this work. Red lines with arrow heads denote 
the excitatory projections mediated by glutamate receptors. Blue solid and dashed lines with round heads represent the inhibitory projections 
mediated by GABAa and GABAb receptors, respectively. Compared with the traditional model of corticothalamic system on absence seizures, the 
basal ganglia are also included in our biophysical model. 
doi:1 0.1 371/journal.pcbi.1 003495.g001 



where D^p is a differential operator representing the dendritic 
filtering of incoming signals, a and fi are the decay and rise times 
of cell-body response to incoming signals, respectively, v^/, is the 
coupling strength between neural populations of type a and type 
b. (I>i,il,t) is the incoming pulse rate from the neural population 
of type b to type a. For simphcity, we do not consider the 
transmission delay among most neural populations in the present 
work. However, since the GABAb functions via second messenger 
processes, a delay parameter T is introduced to its incoming 
pulse rate (i.e., (^^(r,? — t)) to mimic its slow synaptic kinetics. This 
results in a delay differential equation in the fmal mathematical 
description of the BGCT model. Note that the similar modelling 
method has also been used in several previous studies [13,37]. 

In our system, each neural population gives rise to a field 
of pulses, which travels to other neural population at a mean 
conduction velocity v^. In the continuum limit, this type of 
propagation can be well-approximated by a damped wave 
equation [9,33-35,38]: 



(4) 



Here V is the Laplacian operator (the second spatial derivative), 
Va is the characteristic range of axons of type a, and ya = Va/ru 
governs the temporal damping rate of pulses. In our model, only 
the axons of cortical excitatory pyramidal neurons are assumed to 
be sufficiently long to yield significant propagation effect. For other 
neural populations, their axons are too short to support wave 



propagation on the relevant scales. This gives (j)^ = F{ Vc) 
[c = i,r,s,d\,d2,p\,p2,C)- Moreover, as one of typical generalized 
seizures, the dynamical activities of absence seizures are believed 
to occur simultaneously throughout the brain. A reasonable 
simplification is therefore to assume that the spatial activities are 
uniform in our model, which has been shown as the least stable 
mode in models of this class [33,34,36]. To this end, we ignore 
the spatial derivative and set V^ = 0 in Eq. (4). Accordingly, the 
propagation effect of cortical excitatory axonal field (jig is finally 
given by [33,34,36]: 



. d 2 



Kit) = QAt), 



(5) 



where ye = Veli'e- For the population of cortical inhibitory 
interneurons, the BGCT model can be further reduced by using 
Vj = and Q, = Qe, which is based on the assumption that 
intracortical connectivities are proportional to the numbers of 
synapses involved [9,13,33-36]. It has been demonstrated that 
by making these above reductions, the developed BGCT model 
becomes computationally more tractable without significant 
deteriorating the precision of numerical results. 

We then rewrite above equations in the first-order form for 
all neural populations. Following above assumptions, we use Eqs. 
( 1 )-(3) and (5) for modelling the dynamics of excitatory pyramidal 
neurons, and Eqs. (l)-(3) for modelling the dynamics of other 
neural populations. This yields the final mathematical description 
of the BGCT model given as follows: 
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dm 

dt 



dt 



[-^^{t)+F(V,m-2yMt), 



dX(t) 
dt 



-m. 



dX(t) 
dt 



-CYit)-W{t), 



(6) 



(7) 



(8) 



(9) 



where 



Xit) = [ Kit), Frfj it), Vd^it), Vp^ it), Vp^it), V^it), Vrit), Kit) ] ' 
C = oi^, 
Wit) = ia + P)Xit). 



(10) 
(11) 

(12) 



In Eq. (10), the superscript T denotes transposition. The detailed 
expression of Yit) for different neural populations is represented 
by C\, Yiit) and I^2(0> given by 



Yit) = CiY,it)-Y2it), 



(13) 



with 



Ci = 
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1-1(0= 



[^,,F(Ve)MV,0'Piyd2)'nVp0MVp^mV0AF(Vr),F(Vr(t-x))),F(V,)] 



(15) 



Y2it) = 

[ Veil), Vd^ it), Vd2 it), Vp^ it), Vp^it), n(0, Vrit), Kit) - &n ] 



(16) 



Here the variable t in Eq. (15) denotes the GABAb delay and the 
parameter in Eq. (16) represents the constant nonspecific 
subthalamic input onto SRN. 

The parameters used in our BGCT model are compatible with 
physiological experiments and their values are adapted from 
previous studies [9,11,13,30,31,36]. Unless otherwise noted, we 
use the default parameter values listed in Table 1 for numerical 
simulations. Most of the default values of these parameters given in 
Table 1 are based on either their nominal values or parameter 
ranges reported in above literature. A small number of parameters 
associated with the basal ganglia (i.e., Vdjdi> '^pipi ^^'^ ^pid ^'''^ 
adjusted slighdy, but still within their normal physiological 
ranges, to ensure our developed model can generate the stable 
2-4 Hz SWDs under certain conditions. Note that due to lack of 
quantitative data, the coupling strength of the SNr-TRN pathway 
needs to be estimated. Considering that the SNr sends GABAergic 
projections both to SRN and TRN and also both of these two 
nuclei are involved in thalamus, it is reasonable to infer that 
the coupling strengths of these two pathways are comparable. 
For simplicity, here we chose Vrp, = Vj^,, = — 0.035 mV s by default. 
In the following studies, we also change (decrease or increase) 
the value of Vr^,, several folds by employing a scale factor K 
(see below) to examine how the inhibition from the SNr-TRN 
pathway regulates absence seizures. Additionally, during this 
study, several other critical parameters (i.e., v^r, t and Vpjj) are also 
varied within certain ranges to obtain different dynamical states 
and investigate their possible effects on the modulation of absence 
seizures. 

Data analysis 

In the present study, several data analysis methods are 
employed to quantitatively e\'aluat(^ th(- dynamical states as well 
as the properties of SWDs g(;neratc(l by the model. To reveal 
critical transitions between different dynamical states, we perform 
the bifurcation analysis for several key parameters of the model. 
For one specific parameter, the bifurcation diagram is simply 
obtained by plotting the "stable" local minimum and maximum 
values of cortical excitatory axonal fields (i.e., 1^,,) over changes in 
this parameter [11,39]. To this end, all simulations are executed 
for sufficiently long time ( 1 0 seconds of simulation time, after the 
system reaches its stable time series), and only the local minimum 
and maximum values obtained from the latter stable time series 
are used. Using the above bifurcation analysis, we can also easily 
distinguish different dynamical states for combined parameters. 
Such analysis technique allows us to further identify different 
dynamical state regions in the two-parameter space (for example, 
see Fig. 2D). On the other hand, the power spectral analysis is used 
to estimate the dominant frequency of neural oscillations. To do 
this, the power spectral density is obtained from the time series tj)^ 
(over a period of 10 seconds) by using the fast Fourier transform. 
Then, the maximum peak frequency is defined as the dominant 
frequency of neural oscillations. It should be noted that, by 
combining the results of both the state and frequency analysis, we 
can outline the SWD oscillation region that falls into the 2-4 Hz 
frequency range in the two-parameter space (for example, see the 
asterisk region in Fig. 2E). Moreover, we calculate the mean firing 
rates (MFRs) for several key neural populations in some figures. To 
compute the MFRs, all corresponding simulations are performed 
up to 25 seconds and the data from 5 to 25 seconds are used for 
statistical analysis. To obtain convincing results, we carry out 20 
independent simulations with different random seeds for each 
experimental setting, and report the averaged result as the final 
result. Finally, in some cases, we also compute the low and high 
triggering mean firing rates (TMFRs) for SNr neurons. In the 
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Table 1. Default parameter values. 




Symbol 


Description 


Value 


References 




Cortical maximum firing rate 


250 Hz 


[9,11,13] 




Striatum maximum firing rate 


65 Hz 


[30,31] 


Q';r 


SNr maximum firing rate 


250 Hz 


[30,31] 


a;;r 


GPe maximum firing rate 


300 Hz 


[30,31] 


or 


STN maximum firing rate 


500 Hz 


[30,31] 




SRN maximum firing rate 


250 Hz 


[9,11,13] 




TRN maximum firing rate 


250 Hz 


[9,11,13] 


e,, fl,. 


Mean firing threshold of cortical populations 


15mV 


[9,11,13,36] 




IVlean firing threshold of striatum 


19mV 


[30,31] 




IVlean firing threshold of SNr 


lOmV 


[30,31] 




Mean firing threshold of GPe 


9mV 


[30,31] 




Mean firing threshold of STN 


10 mV 


[30,31] 


e. 


Mean firing threshold of SRN 


15mV 


[9,11,13,36] 


0, 


Mean firing threshold of TRN 


15 mV 


[9,11,13,36] 


fee 


Self-coupling strength of EPN 


ImVs 


[11,13] 


— Vei 


Coupling strength from UN to EPN 


1.8mVs 


[11,13] 


Vre 


Coupling strength from EPN to TRN 


0.05 mVs 


[9,13] 




Coupling strength from SRN to TRN 


0.5 mVs 


[9,13] 




Coupling strength from TRN to SRN 


0.4-2mVs 


[9,36] 




Coupling strength from EPN to striatal D1 neurons 


ImVs 


[30,31] 


— 1',;,,/, 


Self-coupling strength of striatal D1 neurons 


0.2 mVs 


[30,31] 




Coupling strength from SRN to striatal D1 neurons 


0.1 mVs 


[30,31] 




Coupling strength from EPN to striatal D2 neurons 


0.7 mVs 


[30,31] 




Self-coupling strength of striatal D2 neurons 


0.3 mVs 


[30,31] 


Vis 


Coupling strength from SRN to striatal D2 neurons 


0.05 mVs 


[30,31] 




Coupling strength from striatal D1 neurons to SNr 


0.1 mVs 


[30,31] 




Coupling strength from GPe to SNr 


0.03 mVs 


[30,31] 




Coupling strength from STN to SNr 


0-0.6mVs 


[30,31] 


- ^pnh 


Coupling strength from striatal D2 neurons to GPe 


0.3 mVs 


[30,31] 


^ ^'P2P2 


Self-coupling strength of GPe 


0.075 mVs 


[30,31] 




Coupling strength from STN to GPe 


0.45 mVs 


[30,31] 


-nP2 


Coupling strength from GPe to STN 


0.04 mVs 


[30,31] 




Coupling strength from SRN to EPN 


1.8mVs 


[9,13] 


V,e 


Coupling strength from EPN to SRN 


2.2 mVs 


[9] 


f£e 


Coupling strength from EPN to STN 


0.1 mVs 


[30,31] 




Coupling strength from SNr to SRN 


0.035 mVs 


[30,31] 




Coupling strength from SNr to TRN 


0.035 mVs 


Estimated 




Cortical damping rate 


100 Hz 


[9,11,13] 


T 


Time delay due to slow synaptic kinetics of GABAg 


50 ms 


[13] 


a 


Synaptodendritic decay time constant 


50s^' 


[9,11,13,36] 


/J 


Synaptodendritic rise time constant 


200s-' 


[9,11,13,36] 


(T 


Threshold variability of firing rate 


6mV 


[9,11,36] 


!*„ 


Nonspecific subthalamic input onto SRN 


2mVs 


[9,11,36] 


Model parameters 
doi:10.1371/journa 


employed in the present study. Unless otherwise noted, we use these default 
.pcbi.1003495.t001 


parameter values for simulations. 





foUowing simulations, we find that the mean firing rate of SNr activation level of SNr in our work (see the Results section). Based 
neurons is increased with the growth of the excitatory coupling on this property, the low and high TMFRs can be determined by 
strength Vy,,^, which serves as a control parameter to modulate the the mean firing rates of SNr neurons occurring at the boundaries 
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Figure 2. Absence seizure activities induced by the slow icinetics of GABAb receptors In TRN. A, B: Bifurcation diagrams of (j)^ as a 
function of the TRN-SRN inhibitory coupling strength — f,,. (A) and the delay parameter t (B), respectively. Four different dynamical states can be 
observed from the time series of (j)^, which are: the saturation state (I), the SWD oscillation state (II), the simple oscillation state (III) and the low firing 
state (IV). C: Typical time series of (j)^, correspond to the above four dynamical states. Here we set T = 50ms and chose i',,. = — 0.48mVs (1), 
Vsr= — 1 mVs (II), v,r = — 1 .48 mV s (111), v„. = — 1.6mVs (IV), respectively. The colors in bifurcation diagrams (A) and (B) correspond to the typical time 
series plotted in (C). D, E: The state analysis (D) and frequency analysis (E) in the ( — v„.,t) panel. Different colors in (D) represent different dynamical 
state regions, corresponding to those dynamical states given in (A), (B) and (C). The yellow asterisk regions surrounded by black dashed lines in (D) 
and (E) represent the SWD oscillation regions falling into the 2-4 Hz frequency range. The other symbols in (D) and (E) are linked to parameter values 
used for different typical time series in (C): I (white filled circle), II (red filled square). III (yellow filled triangle), and IV (purple filled star). For all 
simulations, we set Vpj^=0.3mVs. 
doi:1 0.1 371 /journal.pcbi.1 003495.g002 



of the typical region of 2-4 Hz SWD.S (for example, see the black 
dashed lines in Fig. 3B). 

Numerical simulation 

All network simulations are written and performed under the 
MATLAB environment. The aforementioned dynamical equa- 
tions are integrated by using the standard fourth-order Runge- 
Kutta method, with a fixed temporal resolution of A = 0.05 ms 
[40]. In additional simulations, it turns out that the chosen 
integration step is sufficiendy small to ensure the numerical 
accuracy of our developed BGCT model. The computer 
codes used in the present study are available on ModelDB 



(https:/ /senselab.med.yale.edu/ModelDB/ showmodel.asppmodel = 
152113). The fundamental implementation of the BGCT model 
is provided as supplementary information to this paper (we 
also provide a XPPAUT code for comparison [41]; see Text SI 
and S2). 

Results 

Slow l<inetics of GABAb receptors in TRN create absence 
seizure activities for the BGCT model 

Previous studies have suggested that the slow kinetics of 
GABAb receptors in TRN are a candidate pathological factor 
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Figure 3. Control of absence seizures by the isolated SNr-TRN pathway. A, B: The state analysis (A) and frequency analysis (B) in the 
( — v,,r,v'p,() panel. Here — 1'„- is the inhibitory coupling strength of the TRN-SRN pathway, whereas v^,^- is the excitatory coupling strength of the STN- 
SNr pathway. Different colors in (A) represent different dynamical state regions: the saturation region (I), the SWD oscillation region (II), the simple 
oscillation region (111) and the low firing region (IV). The suppression of SWDs appears to the right of the white dashed line in (A), where the white 
down arrow indicates that the SWD oscillation can be inhibited by decreasing Vp^-. The yellow asterisk region surrounded by black dashed lines in (B) 
denotes the typical 2-4 Hz SWD oscillation region. C, D: Bifurcation diagrams of (j>^. as a function of v^,^- for different — v'.,,. The strengths of the 
inhibitory projections from the TRN to SRN are set as r,.r= — l-20mVs (C) and vv = — 1 .44 mV s (D), respectively. Different colors in (C) and (D) 
represent different dynamical state regions, corresponding to those in phase diagram (A). E: The mean firing rates (MFRs) of several key neural 
populations as a function of r^,^-, with )',., = — 1.44niVs. Here four neural populations are considered: SNr (blue dot), excitatory pyramidal neurons 
(green asterisk), SRN (black circle) and TRN (red square). Two black dashed lines in (E) represent the occurring positions of the low and high triggering 
mean firing rates (TIVlFRs), respectively. F: The low (red filled circle) and high (green filled square) TMFRs as a function of — Vj,. For all simulations, the 
SNr-SRN pathway is artificially blocked (i.e., i'.,,, =OmVs). 
doi:1 0.1 371/journal.pcbi.1 003495.g003 



contributing to the generation of absence seizures both in animal 
experiments and biophysical models of corticothalamic network 
[7,13,42,43]. To explore whether this mechanism also applies 
to the developed BGCT model, we perform one-dimensional 
bifurcation analysis for the inhibitory coupling strength — v^r and 
the delay parameter T, respectively. The corresponding bifurcation 
diagrams and typical time series of (j)^ are depicted in Figs. 2A-2C, 
which reveal that different dynamical sates emerge in our system 



for different values of — Vj,- and T. When the coupling strength 
— Vjr is too weak, the inhibition from TRN cannot effectively 
suppress the firing of SRN. In this case, due to the strong 
excitation from pyramidal neurons, the firing of SRN rapidly 
reaches a high level after the beginning of the simulation. Such 
high activation level of SRN in turn drives the firing of cortical 
neurons to their saturation states within one or two oscillation 
periods (region I). As the couphng strength — Vj,. grows, the 
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inhibition from TRN starts to affect the firing of SRN. For 
sufficiently long T, this causes our model to successively undergo 
two different oscillation patterns. The first one is the SWD 
oscillation pattern, in which multiple pairs of maximum and 
minimum values are found within each periodic complex (region 
II). Note that this oscillation pattern has been extensively observed 
on the EEG recordings of real patients during absence seizures [1]. 
The other one is the simple oscillation pattern, in which only one 
pair of maximum and minimum values appears within each 
periodic: complex (region III). However, if the coupling strength 
— is too strong, the firing of SRN is almost completely inhibited 
by TRN. In this situation, the model is kicked into the low firing 
region and no oscillation behavior can be observed anymore 
(region TV). Additionally, we also find that the model dynamics are 
significanfly influenced by the GABAb delay, and only sufficientiy 
long T can ensure the generation of SWDs in the developed model 
(see Fig. 2B). 

To check whether our results can be generalized within a 
certain range of parameters, we further carry out the two- 
dimensional state analysis in the ( — Vsr,T) panel. As shown in 
Fig. 2D, the whole (— Vv,.,t) panel is divided into four state regions, 
corresponding to those regions identified above. Unsurprisingly, 
we find that the BGCT model can generate the SWD oscillation 
pattern only for appropriately intermediate — Vjr and sufficientiy 
long T. This observation is in consistent with our above finding, 
demonstrating the generafizability of our above results. To 
estimate the frequency characteristics of different oscillation 
patterns, we compute the dominant frequency basc-d on the 
spectral analysis in the (— Vj^jT) panel. For both the simple and 
SWD oscillation patterns, the dominant frequency is influenced by 
T and —Vsr, and increasing their values can both reduce the 
dominant frequency of neural oscillations (Fig. 2E). However, 
compared to — v^,-, our results indicate that the GABAb delay 
may have a more significant effect on the dominant oscillation 
frequency (Fig. 2E). By combining the results in Figs. 2D and 
2E, we roughly outiine the SWD oscillation region that falls into 
the 2-4 Hz frequency range (asterisk region). It is found that most 
of, but not all, the SWD oscillation region is contained in this 
specific region. Here we emphasize the importance of this specific 
region, because the SWDs within this typical frequency range is 
commonly observed during the paroxysm of absence epilepsy in 
human patients [1,2]. 

Why can the slow kinetics of GABAb receptors in TRN induce 
absence seizure activities? Anatomically, the SRN neurons receive 
the TRN signals from the inhibitory pathway mediated by both 
GABAa and GABAb receptors. Under suitable condition, the 
double suppression caused by these two types of GABA receptors 
occurring at different time instants may provide an effective 
mechanism to create multiple firing peaks for the SRN neurons 
(see below). Such firing pattern of SRN in turn impacts the 
dynamics of cortical neurons, thus leading to the generation of 
SWDs. It should be noted that, during the above processes, both T 
and — Vsr play critical roles. In each oscillation period, after the 
GABAA-induced inhibition starts to suppress the firing of SRN 
neurons, these neurons need a certain recovery time to restore 
their mean firing rate to the rising state. Theoretically, if this 
recovery time is shorter than the GABAb delay, another firing 
peak can be introduced to SRN neurons due to the latter 
GABAB-induced inhibition. The above analysis implies that our 
model requires a sufficient long GABAb delay to ensure the 
occurrence of SWDs. However, as described above, too long t is 
also a potential factor which may push the dominant frequency of 
SWDs beyond the typical frequency range. For a stronger — Vsr^ 
the inhibition caused by GABAa is also strong. In this situation. 



it is obvious that the SRN neurons need a longer time to restore 
their firing rate. As a consequent, a relatively longer t is required 
for the BGCT model to ensure the occurrence of SWDs for 

stronger — v,,- (see Fig. 2D). 

These findings provide consistent evidence that our developed 
BGCT model can replicate the typical absence seizure activities 
utilizing previously verified pathological mechanism. Because 
we do not change the normal parameter values for basal gangha 
during above studies, our results may also indicate that, even 
though the basal ganglia operate in the normal state, the abnormal 
alteration within the corticothalamic system may also trigger 
the onset of absence epilepsy. Throughout the following studies, 
we set T = 50ms for all simulations. For this choice, the delay 
parameter T is within the physiological range and modest, allowing 
the generation of SWD oscillation pattern while preserving its 
dominant frequency around 3 Hz in most considered parameter 
regions. It should be noted that, in additional simulations, we have 
shown that by slightiy tuning the values of several parameters our 
developed BGCT model is also powerflil to reproduce many other 
topical patterns of time series, such as the alpha and beta rhythms 
(see Fig. SI), which to a certain extent can be comparable with real 
physiological EEG signals [9,36]. 

Control of absence seizures by the isolated SNr-TRN 
pathway 

Using the developed BGCT model, we now investigate the 
possible roles of basal ganglia in controlling absence seizure 
activities. Here we mainly concentrate on how the activation level 
of SNr influence the dynamics generated by the model. This is 
because, on the one hand, the SNr is one of chief output nucleus of 
the basal ganglia to thalamus, and on the other hand, its firing 
activity has been found to be highly associated with the regulation 
of absence seizures [21,22]. To this end, the excitatory coupling 
strength v^jj is employed to control the activation level of SNr and 
a three-step strategy is pursued in the present work. In this and 
next subsections, we assess the individual roles of two different 
pathways emitted from SNr to thalamus (i.e., the SNr-TRN and 
SNr-SRN pathways) in the control of absence seizures and discuss 
their corresponding biophysical mechanisms, respectively. In the 
final two subsections, we further analyze the combination effects of 
these two pathways on absence seizure control and extend our 
results to more general cases. 

To explore the individual role of the SNr-TRN pathway, we 
estimate both the state regions and frequency characteristics in the 
(— v„.,Vp,j) panel. Note that during these investigations the SNr- 
SRN pathway is artificially blocked (i.e., v,^, =OmVs). With this 
"naive" method, the modulation of absence seizure activities by 
the SNr-SRN pathway is removed and the effect caused by the 
SNr-TRN pathway is theoretically amplified to the extreme. 
Similar to previous results, we find that the whole { — Vsr,Vpji;) panel 
can be also divided into four different regions (Fig. 3A). These 
regions are the same as those defined above. For weak inhibitory 
coupling strength — v^r, increasing the excitatory coupling strength 
Vpjf moves the model dynamics from the SWD oscillation state to 
the saturation state. Here we have to notice that the saturation 
state is a non-physiological brain state even though it does not 
belong to typical seizure activities. In strong —v^r region, the 
suppression of SWDs is observed by decreasing the excitatory 
coupling strength Vpjj, suggesting that inactivation of SNr neurons 
may result in seizure termination through the SNr-TRN pathway 
(Fig. 3 A, right side). For strong enough inhibitory coupling 
strength —Vsr, such suppression effect is rather remarkable that 
sufficientiy low activation of SNr can even kick the network 
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dynamics into the low firing region (compare the results in Figs. 3C 
and 3D). 

The SNr-TRN pathway induced SWD suppression is compli- 
cated and its biophysical mechanism is presumably due to 
competition-induced collision. On the one side, the decrease of 
excitatory coupling strength v^jj inactivates the SNr (Fig. 3E, top 
panel), which should potentially enhance the firing of TRN neurons. 
On the other side, however, increasing the activation level of TRN 
tends to suppress the firing of SRN, which significantly reduces the 
firing of cortical neurons and in turn inactivates the TRN neurons. 
Furthermore, the inactivation of cortical neurons also tends to 
reduce the firing level of TRN neurons. As the excitatory coupling 
strength Vp,( is decreased, the collision caused by such complicated 
competition and information interactions finally leads to the 
inactivation for all the TRN, SRN, and cortical neurons (Fig. 3E, 
bottom panel), which potentially provides an effective mechanism to 
destabilize the original pathological balance within the corticotha- 
lamic system, thus causing the suppression of SWDs. 

Indeed, we find that not only the dynamical state but also the 
oscillation frequency is greatly impacted by the activation level of 
SNr, through the SNr-TRN pathway. For both the simple and 
SWD oscillation patterns, increasing the excitatory strength Vpjj 
can enhance their dominant frequencies. The combined results of 
Figs. 3A and 3B reveal that, for a fixed — Vsr, whether the model 
can generate the SWDs within the typical 2-4 Hz is determined 
by at least one and often two critical values of V'y,,f (Fig. 3B, asterisk 
region). Because the activation level of SNr is increased with the 
growth of Vp,{, this finding further indicates that, due to effect of 
the SNr-TRN pathway, the model might exist the corresponding 
low and high triggering mean firing rates (TMFRs) for SNr 
neurons (Fig. 3E, dashed lines). If the long-term mean firing rate of 
SNr neurons fails into the region between these two TMFRs, the 
model can highly generate typical 2-4 Hz SWDs as those observed 
on the EEG recordings of absence epileptic patients. In Fig. 3F, 
we plot both the low and high TMFRs as a function of the 
inhibitory coupling strength — Vj^. With the increasing of — v,,., the 
high TMFR grows rapidly at first and then reaches a plateau region, 
whereas the low TMFR almost linearly increases during this 
process. Consequentiy, it can be seen that these two critical TMFRs 
approach each other as the inhibitory coupling strength —Vsr is 
increased until they almost reach an identical value (Fig. 3F). 

The above findings indicate that the SNr-TRN pathway 
may play a vital role in controlling the absence seizures and 
appropriately reducing the activation level of SNr neurons can 
suppress the typical 2-A Hz SWDs. The similar antiepUeptic effect 
induced by inactivating the SNr has been widely reported in 
previous electrophysiological experiments based on both genetic 
absence epilepsy rats and tottering mice [21-23,26]. Note that, 
however, in literature such antiepUeptic effect by reducing the 
activation of SNr is presumed to be accomplished through the 
indirect SNr-TRN pathway relaying at superior coUiculus [21,22]. 
Our computational results firstly suggest that such antiepileptic 
process can be also triggered l)y the direct SNr-TRN GABAergic 
projections. Combining these results, we postulate that for real 
absence epileptic patients both of these two pathways might 
work synergistically and together provide a stable mechanism to 
terminate the onset of absence epilepsy. 

Control of absence seizures by the isolated SNr-SRN 
pathway 

We next turn on the SNr-SRN pathway and investigate whether 
this pathway is also effective in the control of absence seizures. 
Similar to the previous method, we artificially block the SNr-TRN 



pathway (i.e., Vrp, =OmVs) to enlarge the effect of the SNr-SRN 
pathway to the extreme. Fig. 4A shows the two-dimensional state 
analysis in the (— v„-,v,,|r) panel, and again the whole panel is 
divided into four different state regions. Compared to the results in 
Fig. 3A, the suppression of SWDs appears in a relatively weaker 

— Vsr region by increasing the excitatory coupling strength v^, j . This 
finding suggests that the increase in the activation of SNr can also 
terminate the SWDs, but through the SNr-SRN pathway. For 
relatively weak — tv within the suppression region, the SNr-SRN 
pathway induced suppression is somewhat strong. In this case, the 
high activation level of SNr directly kicks the network dynamics into 
the low firing region, without undergoing the simple oscillation state 
(Fig. 4C2 and compare with Fig. 4C3). Note that this type of state 
transition is a novel one which has not been observed in the SWD 
suppression caused by the SNr-TRN pathway. For relatively strong 

— v„. within the suppression region, the double peak characteristic 
of SWDs generated by our model is weak. In this situation, as the 
inhibitory coupling strength — Vsr is increased, we observe that the 
network dynamics firstiy transit from the SWD oscillation state to 
the simple oscillation state, and then to the low firing state (Fig. 4C3). 

To understand how the SNr-SRN pathway induced SWD 
suppression arises, we present the mean firing rates of several key 
neural populations within the corticothalamic system, as shown 
in Fig. 4D. It can be seen that increasing the strength Vpjj 
significanfly improves the activation level of SNr (Fig. 4D, 
top panel), which in turn reduces the firing of SRN neurons 
(Fig. 4D, bottom panel). The inactivation of SRN neurons further 
suppresses the mean firing rates for both cortical and TRN 
neurons (Fig. 4D, bottom panel). These chain reactions lead to the 
overall inhibition of firing activities in the corticothalamic system, 
which weakens the double peak shaping effect due to the slow 
kinetics of GABAb receptors in TRN. For strong Vp,j, such 
weakening effect is considerable, thus causing the suppression of 
SWDs. Our results provide the computational evidence that high 
activation of SNr can also effectively terminate absence seizure 
activities by the strong inhibition effect from the SNr-SRN pathway. 
Compared to the SWD suppression induced by the SNr-TRN 
pathway, it is obvious that the corresponding biophysical mecha- 
nism caused by the SNr-SRN pathway is simpler and more direct. 

Moreover, our two-dimensional frequency analysis indicates 
that the dominant frequency of neural oscillations depends on the 
excitatory coupling strength Vpjj (see Fig. 4B). For a constant — Vsr, 
progressive increase of Vpji; reduces the dominant frequency, 
but not in a very significant fashion. Thus, we find that almost all 
the SWD oscillation region identified in Fig. 4A falls into the 
typical 2-4 Hz frequency range (Fig. 4B, asterisk region). Unlike 
the corresponding results presented in previous subsection, the 
combination results of Figs. 4A and 4B demonstrate that the 
BGCT model modulated by the isolated SNr-SRN pathway only 
exhibits one TMFR for SNr neurons. For a suitably fixed strength 

— Vsr, the generation of SWDs can be highly triggered when the 
mean firing rate of SNr neurons is lower than this critical firing 
rate (Fig. 4D, dashed line). With the increasing of —Vsr, we 
observe that this TMFR rapidly reduces from a high value to a low 
value (Fig. 4E). Note that this decreasing tendency is in contrast 
with our previous finding based on the model modulated by the 
isolat(xl SNr-TRN patln\ay (compared to the results in Fig. 3F). 

Takc'n together, these observations suggest that increasing the 
acti\ ation level of SNr neurons also significantly influences the 
dynamics of the corticothalamic system and causes the suppression 
of absence seizure activities. To the best of our knowledge, this is a 
new finding that underscores the importance of the direct 
inhibitory SNr-SRN pathway in controlling and modulating 
absence seizure activities. It is reasonable to believe that several 
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Figure 4. Control of absence seizures by the isolated SNr-SRN pathway. A, B: The state analysis (A) and frequency analysis (B) in the 
( — I'll ,V;)iC) panel. Here — iv is the inhibitory coupling strength of the TRN-SRN pathway, whereas v^,^- is the excitatory coupling strength of the STN- 
SNr pathway. Similar to previous results, four different dynamical state regions are observed: the saturation region (I), the SWD oscillation region (II), 
the simple oscillation region (III) and the low firing region (IV), which correspond to those defined in Fig. 2 (D). The region between two white dashed 
lines in (A) represents the suppression region of SWDs, where the white up arrow indicates that the SWD oscillation can be inhibited by increasing 
Vp,;. The yellow asterisk region surrounded by black dashed lines in (B) denotes the typical 2-4 Hz SWD oscillation region. C: Two different types of 
SWD suppressions caused by linearly increasing Vp,{. Top (C,): The value of 1'^,^- as a function of time. Middle (C2): Corresponding (j)^ trace for 
Vj,. = — 0.72mVs. Bottom (C3): Corresponding (ji^, trace for v„= — 0.92mVs. D: The MFRs of several key neural populations as a function of 1',,,^, with 
IV = — 0.72mVs. Here four neural populations are considered: SNr (blue dot), excitatory pyramidal neurons (green asterisk), SRN (black circle) and 
TRN (red square). The black dashed line in (D) represents the occurring position of TMFR. E: The TMFR as a function of — v.,, . Note that the SNr-TRN 
pathway is artificially blocked (i.e., iv^, =OmVs) for all simulations. 
doi:1 0.1 371 /journal.pcbi.1 003495.g004 



other external factors, which are able to enhance the activation 
level of SNr, may also result in the termination of absence seizures 
due to the similar mechanism. 

Competition-induced bidirectional control of absence 
seizures by the basal ganglia 

So far, we have confirmed that the absence seizure activities 
generated by the BGCT model can be inhibited by either the 



isolated SNr-TRN pathway or the isolated SNr-SRN pathway, 
through different biophysical mechanisms. In real brain, however, 
both of these pathways should be available and work together at 
the same time. Thus, an important and naturally arising question 
is whether these two types of seizure control can coexist in the 
same network, and if possible, whether this feature can be 
maintained in certain range of parameters. To address this issue, 
we introduce a scale factor K and set the coupling strength 
v,y,, = KVy,^ . By utilizing this method, we can flexibly control the 
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Figure 5. Control of absence seizures by the combination effects of the SNr-TRN and SNr-SRN pathways. A, B: The state analysis (A) and 
frequency analysis (B) in the (K,Vp^^) panel. Here K is the scale factor, and v^,^ is the excitatory coupling strength of the STN-SNr pathway. Unlike 
previous results, only three dynamical state regions are observed in the phase diagram (A): the SWD oscillation region (II), the simple oscillation 
region (III) and the low firing region (IV), corresponding to the dynamical states defined in Fig. 2(D). For relatively weaker scale factor K, both increase 
and decrease in the activation level of SNr can inhibit the SWDs (double arrow, bidirectional suppression); whereas for sufficiently strong K, only 
reducing the activation level of SNr suppresses the SWDs (single arrow, unidirectional suppression). In (A), the white dashed line represents the 
boundary of suppression region, and the red dashed line stands for the demarcation between the bidirectional and unidirectional suppression 
regions. The yellow asterisk region surrounded by black dashed lines in (B) denotes the SWD oscillation region that falls into the 2-4 Hz frequency 
range. C: The MFRs of several key neural populations as a function of Vp,^, with the scale factor K = 0.6. Here four neural populations are considered: 
SNr (blue dot), excitatory pyramidal neurons (green asterisk), SRN (black circle) and TRN (red square). Two black dashed lines in (C) represent the 
occurring positions of low and high triggering mean firing rates (TMFRs), respectively. D: The low (red filled circle) and high (green filled square) 
TMFRs as a function of K. For all simulations, the coupling strength of the TRN-SRN pathway is set as iv = — 1-08 mVs. 
doi:1 0.1 371 /journal.pcbi.1 003495.g005 



relative coupling strength between the SNr-TRN and SNr-SRN 
pathways and discuss the combination roles of these two inhibitory 
pathways in detail. 

Fig. 5 A shows the state analysis in the {K,Vp^;;) panel with 
Vsr = — 1 .08 mV S. Unlike the previous results, here we only 
discover three dynamical state regions, which correspond to: the 
SWD oscillation state (II), the simple oscillation state (III), and the 
low firing state (IV). The disappearance of the saturation state is at 
least due to the following two reasons: (1) the double suppression 
from the SNr-TRN and SNr-SRN pathways and (2) the relatively 
strong inhibitory effect from TRN to SRN. As shown in Fig. 5A, 
the phenomenon of the SWD suppression appears in the strong 
K region. For relatively weak K within this suppression region, 
both increasing and decreasing the activation of SNr neurons from 
the normal level effectively inhibit the generation of SWDs (see 
Fig. 5A). Such bidirectional suppression behavior can be attributed 
to the effective competition between the SNr-TRN and SNr-SRN 
pathways. As the scale factor K is increased, the enhancement of 
the SNr-TRN pathway breaks the original competition balance 



between these two inhibitory pathways. During this process, the 
inhibition from the SNr-TRN pathway progressively dominates 
the model dynamics. Accordingly, for sufficiently strong K, the 
suppression of SWDs is only found by lowing the activation of SNr 
(Fig. 5A), which is consistent with our previous critical observation 
given in Fig. 3A. 

However, as described above, the increase in the excitatory 
coupling strength v^,f significantiy reduces the dominant frequen- 
cy of SWDs, through the SNr-TRN pathway. Therefore, although 
enhancing the activation level of SNr neurons cannot inhibit the 
SWDs direcdy, it tends to push the dominant frequency of SWDs 
below 2 Hz (Fig. 5B). By combining the results of Figs. 5A and 5B, 
we successfully outline the SWD oscillation region that falls into 
the 2-4 Hz frequency range (asterisk region). In this typical SWD 
region, the model exhibits both the low and high TMFRs of the 
SNr neurons for a fixed K (Figs. 5B and 5C). The generation of 
SWDs within the typical frequency range can be highly triggered if 
the mean firing rate of SNr neurons is between these two critical 
TMFRs (Fig. 5C). Nevertheless, due to the competition between 
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the SNr-TRN and SNr-SRN pathways, the firing activities of 
the TRN, SRN, and cortical neurons become more complicated, 
compared to the cases induced by any isolated pathway. For a 
constant K, such competition creates a bell-shaped MFR curve for 
each key neural population by tuning the excitatory coupKng 
strength v^,^', (see Fig. 5C). To further investigate how the scale 
factor impacts the low and high TMFRs, we plot these two 
TMFRs as a function of isT in Fig. 5D. With the increasing 
of relative strength K, both of these two critical TMFRs are 
rapidly changed but in diiferent fashions. The high TMFR is 
increased from a relatively low value to saturation, whereas the 
low TMFR is reduced from a relatively high value to 0 (see 
Fig. 5D). Obviously, such opposite tendencies of these two TMFRs 
are attributed to the combination effects of the SNr-TRN and 
SNr-SRN pathways. 

Furthermore, we also find that both the suppression of SWDs 
and the typical 2-4 Hz SWD region are shaped by the strength of 
inhibitory projections from the TRN to SRN. In Figs. 6A and 6B, 
we perform a series of two-dimensional state and frequency 
analysis in the {K,Vp^i) panel for different values of — v,,.. When 
the inhibitory coupling strength — v„. is too weak, the SWDs 
generated by the corticothalamic system is mainly controlled by 
the SNr-SRN pathway. In this case, the suppression of SWDs is 
observed in the intermediate K region and only increasing the 
activation level of SNr can effectively inhibit the SWDs (Fig. 6A). 
As the coupling strength — v^r is increased, the inhibition from 



SNr-TRN pathway starts to influence the model dynamics. This 
introduces the competition between the SNr-SRN and SNr-TRN 
pathways, leading to the emergence of the bidirectional suppres- 
sion of SWDs for intermediate K (Fig. 6A). It is obvious that 
the higher the coupling strength — v,,-, the stronger the inhibition 
effect caused by the SNr-TRN pathway. With increasing the value 
of — Vjr, such strengthened inhibition significantiy moves the 
boundary of the low TMFR toward higher values of v^,(, and thus 
notably shrinks the region of SWDs within the typical frequency 
range of 2-4 Hz (Fig. 6B). 

These above observations emphasize the importance of the 
combination role of both the SNr-TRN and SNr-SRN pathways 
on the control of absence seizure activities. Quite remarkably, 
we observe that the bidirectional suppression of SWDs emerges 
under suitable conditions, that is, both increasing and reducing the 
activation levels of SNr neurons effectively suppress the SWDs. 
Such bidirectional suppression is determined and modulated by 
both the relative strength of these two inhibitory pathways and the 
strength of inhibitory projections from the TRN to SRN. This 
novel finding indicates the possible bidirectional control of absence 
seizures by the basal ganglia, which is induced by the competition 
between the SNr-TRN and SNr-SRN pathways. In additional 
simulations, we have found that several other factors, which can 
effectively change the activation level of SNr, can also lead to the 
bidirectional suppression of SWDs due to the similar mechanism, 
further demonstrating the generality of our results. 
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Figure 6. Effect of inhibitory coupling strength — v„ on the control of absence seizures by the SNr-TRN and SNr-SRN pathways. Tow- 
dimensional state analysis (A) and corresponding frequency analysis (B) in the (K,Yp^^) panel for different values of i'„.. Here only three dynamical state 
regions are observed in (A): the SWD oscillation region (II), the simple oscillation region (III) and the low firing region (IV). In (A,)-(A4), the regions 
marked by red diamonds denote the whole suppression regions of SWDs, the white dashed lines represent the boundaries of suppression regions, 
and the red dashed lines stand for the demarcations between the bidirectional (double arrow) and unidirectional (single arrow) suppression regions. 
In (B,)-(B4), the yellow asterisk regions surrounded by black dashed lines denote the typical 2-4 Hz SWD oscillation regions. From left to right, the 
strengths of inhibitory projections from the TRN to SRN are: iv=-lmVs (A,, B,), v,T=-1.04mVs (A2, B2), i'„.= -1.12mVs (A3, B3), and 
v„. = — 1.24 mVs (A4, B4), respectively. For better showing and comparing the bidirectional suppression regions among different subfigures, here we 
mainly consider the scale factor for all values of v.,, within the same and small interval from 0 to 1. As an additional comparison, both the state 
analysis and frequency analysis for v,r= — 1.20mVs in a relatively larger K interval is given in Fig. S2. 
doi:1 0.1 371/journal.pcbi.1 003495.g006 
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Bidirectional control of absence seizures by the basal 
ganglia might be extended to other pathological factors 

In above subsections, we focused on one specific pathological 
factor and have computationally shown that the basal ganglia may 
bidirectionaUy control and modulate the typical absence seizure 
activities (i.e., the SWDs) induced by the slow synaptic kinetics 
of GABAg in TRN. A natural question to ask is whether such 
bidirectional control feature caused by the basal ganglia is a 
generalized regulatory mechanism for absence seizures. We argue 
that this might be true. To check this postulation, at least one 
additional SWD generation mechanism should be introduced into 
our model, and we need to examine whether the bidirectional 
control of absence seizures by the basal ganglia is also available for 
the new pathological factor. In literature, there are several other 
theories that are associated with the generation mechanisms of 
SWDs. A boldly accepted one is related to the transmission delay 
between the cerebral cortex and thalamus and, specifically, it has 
been found that suitably choosing such transmission delay can 
drive the corticothalamic system produce the SWDs [11,33]. To 
apply this pathological factor in our model, here we block the 
GABAb pathway from TRN to SRN, and consider a bidirectional 
?o/2 transmission delay between the cerebral cortex and thalamus, 
as that used in previous studies [11,33]. Additionally, several 
coupling strengths within the corticothalamic loop are also needed 
to be adapted, because such transmission delay induced SWDs 



require strong interactions between the cerebral cortex and 
thalamus [1 1,33]. The new added and modified model parameters 
that we used in this subsection are as follows [11,33]: ro = 80ms, 
V(..5 = 3.2mVs, Vst> = 3.4mVs, v,.t> = 1.6mVs, Vj,.= — 1.76mVs 
and (/>„ = 8.0mVs. For simplicity, we term the current model as 
the modified model in the following studies. 

Figs. 7A and 7B show an example pair of state analysis and 
frequency analysis in the (A^,Vp,f) panel, respectively. As expected, 
due to the competition between the SNr-TRN and SNr-SRN 
pathways, we observe the significant bidirectional control feature 
for intermediate scale factor K (Fig. 7A, the region between 
dashed lines). In this bidirectional region, both enhancing and 
lowing the excitatory coupling strength v^,j push the model 
dynamics from the SWD oscillation state into the simple oscillation 
state (Fig. 7C), thus inhibiting the generation of SWDs. 
This finding supports our above hypothesis that under suitable 
conditions the basal ganglia may control and modulate the 
absence seizure activities bidirectionaUy. However, compared to 
the results in Fig.5, we fmd that the bidirectional region appears in 
a relatively larger K region for the modified model, and increasing 
the coupling strength Vpjj cannot kick the model dynamics into the 
low firing state as well. This is not so surprising because these two 
SWD generation mechanisms that we used are similar but not 
completely identical. As we introduced above, the SWDs induced 
by the corticothalamic loop transmission delay require relatively 
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Figure 7. An example illustrating that the bidirectional control of absence seizures by the basal ganglia is also available for the 
corticothalamic loop transmission delay induced SWDs. A, B: Tow-dimensional state analysis (A) and corresponding frequency analysis (B) in 
the (K,Vp^^) panel. Here K is the scale factor, and v,,,^- represents the excitatory coupling strength of the STN-SNr pathway. Four different dynamical 
states are observed in (A): the saturation state (I), the SWD oscillation state (II), the simple oscillation state (111) and the low firing state (IV), which also 
correspond to previous figures. In (A), the region between two red dashed lines denotes the main bidirectional suppression region of SWDs, where 
the double arrow represents both increasing and decreasing the excitatory coupling strength \',„^ can inhibit the SWDs. In (B), the yellow asterisk 
region surrounded by black dashed lines denotes the typical 2-4 Hz SWD oscillation region. (C) Three typical time series of for different values of 



Vp^-^, with .^;; = 1.3. Here we choose Vp,f = 
fo = 80 ms, v„ = 3.2 mV s, vv = 3.4 mV s, 
doi:1 0.1 371/journal.pcbi.l 003495.g007 



= 0.09mVs (top), v,„{ = 0.6mVs (middle) and v^,; = 
v„ = 1 .6 mV s, v„ = - 1 .76 mV s and (j)„ = 8.0 mV s. 



3.0 mVs (bottom), respectively. In all simulations, we set 
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stxonger interactions between the cerebral cortex and thalamus, 
which essentially weaken the inhibitory effect from the SNr 
neurons. This might lead that the firing of SRN neurons cannot be 
fully suppressed even when the activation level of SNr reaches its 
saturation state. Another important finding that we discover here 
is that the dynamics of the modified model become complicated 
for large K and strong Vpjj (see Fig. 7A, upper right). In additional 
simulations, we further perform a series of state analysis for the 
modified model using the same group of parameter values but 
different random initial conditions (see Fig. S3). The correspond- 
ing results indicate that the modified model shows bistability 
(the simple oscillation state or the saturation state) in the large K 
and strong Vp,f region, and final model dynamics significantiy 
depend on the initial conditions. 

In conclusion, these findings further stress the combination role 
of the inhibitory SNr-TRN and SNr-SRN pathways on the control 
of absence seizure activities. By combining all of our results, we 
postulate that the bidirectional control by the basal ganglia is 
possible a generalized regulatory mechanism for absence seizures 
and may be extendable to other pathological factors, even though 
the detailed bidirectional control behaviors may not be completely 
identical for different pathological factors. 

Discussion 

Using a mean-field macroscopic model that incorporates the 
basal ganglia, cerebral cortex and thalamus, we presented here 
the first investigation on how the basal ganglia control the absence 
epilepsy through the projections direcdy emitted from the SNr 
to several key nuclei of thalamus. Through simulations, we 
demonstrated that the absence seizure activities induced by the 
slow synaptic kinetics of GABAb in TRN can be inhibited 
by either the isolated SNr-TRN pathway or the isolated SNr- 
SRN pathway via different biophysical mechanisms. More 
importandy, our results showed that under certain conditions 
these two types of seizure control can coexist in the same network, 
suggesting that both decreasing and increasing the activation 
levels of SNr may considerably suppress the generation of SWDs. 
Theoretically, such bidirectional control of absence seizures 
by basal ganglia is due to the effective competition between the 
SNr-TRN and SNr-SRN pathways, which might be a generalized 
mechanism for regulating absence seizure activities and can be 
extended to other pathological factors. In addition, our detailed 
frequency analysis also indicated that, depending on difierent 
system conditions, the developed model may exist low, high or 
both TMFRs for the SNr neural population for triggering the 
typical 2-4 Hz SWDs. 

These results are partiy in agreement with former experimental 
observations. Previously, experimental studies based on electro- 
physiological recordings have established the linkage that reducing 
the activation of SNr neurons from the normal level can effertively 
suppress the SWDs in different rodent animal models [21-23,26]. 
Such antiepileptic effect was supposed to be attributed to 
the indirect pathway of the SNr to TRN relaying at superior 
coUiculus. The results presented in this work also demonstrated 
that decreasing the SNr activity is an effective approach that 
terminates the SWDs. However, it is important to note that in our 
model the similar antiepileptic effect is triggered by the direct 
GABAergic jjrojcctions from SNr to TRN. Presumably, this is 
because the SNr has overall inhibitory impacts on TRN via both 
indirect and direct pathways. In the brain of absence epileptic 
patients, both of these two pathways might work together and 
provide a stable and endogenous mechanism to terminate the 
paroxysm of absence epilepsy. 



Our model further makes prediction that increasing the 
activation of SNr from the normal level may also suppress 
SWDs. In previous experimental studies, there still lacks sufficient 
evidence to support this viewpoint. We speculate that this might be 
because the relative strength of the SNr-TRN and SNr-SRN 
pathways for rodent animals is generally high or at least not too 
low. Under this condition, the BGCT system for most rodent 
animals does not operate in the bidirectional control region, thus 
the suppression of SWDs caused by activating SNr is difficult to 
be observed in normal experiments. According to our results, the 
activation of SNr induced SWD suppression might appear by 
suitably tuning the relative strength between these two pathways. 
Further experiments based on animal models will be necessary to 
validate this prediction and characterize the detailed nature of the 
SWD suppression induced by the SNr-SRN pathway. Even so, the 
above prediction from our computational study might provide an 
alternative approach for terminating absence seizures. 

An interesting and important question is: can real brain utilize 
this bidirectional modulation mechanism to control the paroxysm 
of absence epilepsy? Our answer is that it is possible, if the real 
brain has some mechanisms to automatically adjust the balance 
between the SNr-TRN and SNr-SRN pathways. Theoretically, 
there are several possible biological mechanisms and one of them 
is discussed as follows. Experimental data have uncovered that 
synapses conduct signals in an unreliable fashion, which is due 
to the probabilistic neurotransmitter release of synaptic vesicles 
[44-46]. It has been shown that the transmission failure rate at a 
given synapse generally tends to exceed the fraction of successful 
transmission, and in some specific cases it can be even higher 
than 0.9 [44-4r6]. Interestingly, recent studies indicated that such 
synaptic unreliability may play critical functional roles in neural 
processing and computation [47-49]. For patients with absence 
seizures, the suitable competition between the SNr-TRN and SNr- 
SRN pathways can be theoretically achieved by properly tuning 
the synaptic transmission rate of these pathways. This might be an 
important underlying mechanism and has significant functional 
advantages, because it does not require the changes of related 
anatomical structures and connection densities in the brain. 
However, we should notice that such tuning is not easy. 
Specifically, it requires the cooperation among neurons and needs 
to change the synaptic transmission rates collectively for most of 
relevant synapses. From the functional perspective, this mecha- 
nism may be associated with the self-protection ability of brain. 
After a long time of evolution, it is reasonable to suppose that 
our brain might have the abilities to use this t)pe of plasticity- 
like mechanism to achieve complicated self-protection function 
during absence seizures [49]. Nevertheless, additional well- 
designed experiments are still needed to test whether our proposed 
hypothesis is correct. 

The current results highlight the functional roles of basal ganglia 
in the control of absence seizures and might offer physiological 
implications on different aspects. First, the termination of absence 
seizures by activating SNr neurons through the SNr-SRN path- 
way might inspire the treatment of refractory absence seizures. 
Although the absence epilepsy is one typical benign epilepsy, 
a considerable proportion of patients yet may fail to achieve 
freedom from absence seizures and become refractor^' to multiple 
antiepileptic drugs [50,51]. For those patients, one possible reason 
is that the strengths of their direct and indirect SNr-TRN 
pathways are somewhat weak, compared to the other patients. In 
this case, reducing the activation level of SNr may lose efficacy in 
the suppression of SWDs and, contrarUy, increasing the firing of 
SNr neurons may stop the absence seizures. Second, our results 
might be generalized to the GPi neurons. In addition to the SNr, 
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the GPi is also an important output structure of basal ganglia. The 
SNr often works in unison with GPi, since they have closely related 
inputs and outputs, as well as similarities in cytology and function 
[30-32]. It is thus reasonable to infer that the activation level of 
GPi might serve a similar bidirectional role in modulating the 
absence seizures. Third, the results presented in this study may 
provide new insight into the deep brain stimulation therapy on 
absence seizures [52]. In previous studies, it has been demon- 
strated that absence seizures can be inhibited by suitably applying 
the deep brain stimulation to STN. Our results not only support 
the traditional \'i(;wpoint, but also indicate that the SNr may be a 
more direct therapeutic target, compared to STN. We thus infer 
that the inhibition of typical absence seizure activities can also 
be accomplished by appropriately stimulating the SNr neurons. 
Finally, the main goal of this study was to explore the control role 
of basal ganglia in absence seizures, but not limited to this specific 
epilepsy. Several other types of epilepsies, such as the juvenile 
myoclonic epilepsy [53] and the generalized tonic-clonic epilepsy 
[54], are also highly associated with the corticothalamic system 
and mediated by GABA receptors. If our above findings on 
absence seizures capture the real fact, such bidirectional control 
feature by basal ganglia may be also available for these types of 
epilepsies. 

Although our developed model is powerful enough to 
suggest some functional roles of basal ganglia in controlling 
and modulating absence seizures, we have to admit that this 
biophysical model is simplified and idealized. The limitations 
of this model and possible extensions in future studi(;s need to 
be discussed. First, more detailed models based on the spiking 
neurons are able to introduce much more complexities, such as 
ionic dynamics, firing variability and connection property, to 
the system. In previous studies, it has been observed that ion 
concentrations inside the cell and in the extracellular space 
as well as the firing dynamics of neurons are changed during 
epileptic seizures [55-57]. It has been also found recently that the 
connection property of neural systems in difiFerent scales (for 
example, neuronal networks at microscopic scales and brain 
region networks at macroscopic scales) also contribute to the 
generation of seizure-like activities [58,59]. Therefore, we cannot 
exclude other possible interesting results by using more detailed 
modelling methodologies. Moreover, due to lack of necessary data, 
our model does not include the indirect pathway from the SNr to 
TRN, relaying at the intermediate and superficial layers of the 
caudal superior coUiculus [60]. Although we have inferred that 
both the direct and indirect SNr-TRN pathways play the similar 
effects on the control of absence seizures, further studies based on 
detailed anatomical data need to be investigated in the future 
work. Finally, it should be emphasized that in the present study we 
mainly focus on the SWDs induced by the slow kinetics of 
GABAb receptors in TRN. Even though we also showed that the 
similar bidirectional control behavior can b(- obser\ ed for another 
pathological factor (i.e., the corticothalamic loop transmission 
delay), more detailed analysis on this aspect is stiU needed. As 
we know, several other pathological factors, such as the increased 
T-type Ca^"*" current in thalamocortical neurons [2], may also lead 
to the generation of SWDs. Nevertheless, it is still unknown and 
deserves to be further examined whether our proposed bidirec- 
tional control by basal ganglia is also available for these other 
pathological factors. 

Theoretically, a better choice for mimicking the slower 
dynamics of GABAb currents should use relatively smaller values 
of a and p at the relevant synapses. However, we have to notice 
that the mean-filed theory used in the present study is under an 
assumption of uniform a and P parameters for all incoming 



connections of one neural population (see Eq. 3). To our 
knowledge, introducing different a and j6 to different types of 
connections wiU make such simple form of differential operator 

presented in Eq. 3 become invalid, thus making the dendritic 
filtering of incoming signals mathematically intractable in the 
current mean-field framework. Considering this, we employed a 
delay parameter t to represent the order of magnitude of rise and 
decay time of GABAb currents in our model, as the method used 
in a recent neural mean-field modelling study [13]. Essentially, 
more detailed models based on spiking neurons will allow us 
to consider difiFerent dynamical processes for difiFerent types of 
synapses [49,61]. We are currentiy trying to extend our above- 
identified results to a relevant spiking neural network and the 
corresponding results wiU be presented elsewhere. 

It should be also pointed out that traditional modelling studies 
on absence seizures mainly focus on the possible generation 
mechanism of the typical absence seizure activities, i.e., SWDs, 
within the corticothalamic system, but littie research addresses the 
regulatory mechanisms of absence seizures by basal gangUa. 
Inspired by previous experimental observations [21-23,26] and 
recent modelling studies on Parkinson's disease [30,31], here we 
introduced the basal ganglia to the traditional mean-field model of 
corticothalamic system proposed by Robinson and his collabora- 
tors, and firstiy explored how the basal ganglia control the SWDs 
through the computational approach. Theoretically, several other 
models that describe the macro-scale dynamics about the activity 
of neural ensembles, such as the WUson-Cowan model [62-64] 
and the models based on dynamic causal modelling [65-67], 
can be also used to deal with the similar problem but may need to 
re-estimate parameters according to different model assumptions. 

In summary, we have performed mechanistic studies and 
investigated the detailed roles of basal ganglia in the control of 
absence seizures. We have computationally demonstrated that the 
SWDs generated by the developed model can be terminated and 
modulated by both decreasing and increasing the activation levels 
of SNr. Our results provide the first evidence on the bidirectional 
control of absence seizures by basal ganglia. This finding might 
deepen our conventional understanding about the functional roles 
of basal ganglia in the controlling and modulating of absence 
seizures. For patients with absence epilepsy, these results in turn 
indicate that the loss of several basal ganglia functions, especially 
the functions related to the SNr neurons, might further aggravate 
the attacks of absence epilepsy. We hope that predictions from our 
systematic model investigation can not only inspire testable 
hypotheses for future electrophysiological experiments but also 
provide additional therapeutic strategies for absence seizures. 

Supporting Information 

Figure SI Several other typical time series of (j>^ (left 
frames) and their corresponding spectra (right frames) 
generated by our developed BGCT model. To a certain 
extent, these model time series are comparable with real 
physiological EEG signals: eyes-open (A) and (B), alpha rhythm 
(C), beta rhythm (D), coexistence of alpha and beta rhythms (E), 
and polyspike and wave (F). Here we adjust the excitatory 
corticothalamic coupling strength v,^ and delay parameter T for 
reproducing these time series. The detailed parameter values 
used in our simulations are: V5e = 0.6mVs and T = 32ms (A), 
Vj(, = 1.1 mVs and t= 16ms (B), Vje = I.2mVs and t= 16ms (C), 
Vje = 1.25mVs and T = 3ms (D), Vjc = 1.22mVs and T = 16ms 
(E), and = 2.2 mV s and t = 70 ms (F), respectively. In addition, 
we also introduce a certain level of gaussian white noise into 
(j>„, with the mean <((^„ ) = 2.0 mV s and standard deviation 
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(t(0„) = 0.0141 mVs'^^. Note that more types of model time series 
might be observed by further tuning these critical parameter 
values. 
(TIF) 

Figure S2 An example representation of the bidirec- 
tional control of absence seizures by the basal ganglia in 
a relatively larger scale factor interval. The state analysis 
(A) and frequency analysis (B) in the (K,Vp^^) panel, with the 
inhibitory couphng strength Vjr = — 1 .2 mV s. Here K is the scale 
factor, and is the excitatory coupling strength of the STN-SNr 
pathway used to control the activation level of SNr neurons. 
Similar to the results in Figs. 5A and 6A, only three dynamical 
state region.s are ob.served in the phase diagram (A): the SWD 
oscillation region (II), the simple oscillation region (III) and the low 
firing region (IV). In (A), the region marked by red diamond 
denotes the whole suppression regions of SWDs, the white dashed 
line represents the boundary of suppression region, and the red 
dashed line stands for the demarcation between the bidirectional 
(double arrow) and unidirectional (single arr(n\) suppression 
regions. In (B), the red asterisk region surrounded by dashed lines 
denotes the typical 2-4 Hz SWD oscillation region. Compared to 
the results shown in Figs. 5 and 6, here we consider a relatively 
larger scale factor interval from 0 to 2.3. 
(TIF) 

Figure S3 A series of two-dimensional state analysis in 
the (K,Vpji^) panel for the modified model. In (A) (H), we 
use the same group of parameter values but different random 
initial conditions for simulations. Similar to previous results, four 

References 

1. Panayiotopoulos CP (1997) Absence epilepsies. In: Engel J, Pedley TA, eds. 
Epilepsy: a comprehensive textbook. Philadelphia: Lippincott-Raven. pp. 2327— 
2346. 

2. Crunelli V, Leresche N (2002) Childhood absence epilepsy: genes, channels, 
neurons and networks. Nat Rev Neurosci 3:37 1-382. 

3. Marescaux C, Vergnes M (1995) Genetic absence epilepsy in rats from 
Strasbourg (GAERS). ItalJ Neurol Sci 16:113-118. 

4. Cocncn AN4, Van Luijtclaar EL (2003) Genetic animal models for absence 
epilepsy: a review of the WAG/Rij strain of rats. Behav Genel 33:635 65:"). 

5. rimol'eex' 1. Sleriadc M (2004) Xeoeortieal seizures: initiation, development anel 
cessation. Neiirosri 123:299 336. 

6. Lytton WW (2008) (Jompiiler modelling of epilepsy. Nat Rev Neurosei 9:626- 
637. 

7. Destexhe A (1998) Spikc-and-wave oscillations based on the properties of 
GABAb receptors. J Neurosci 18:9099-9111. 

8. Destexhe A (1999) Can GABAa conductances explain the fast oscillation 
frequency of absence seizures in rodents? Eur J Neurosci 11:2175-2181. 

9. Robinson PA, Rennie CJ, Rowe DL (2002) Dynamics of large-scale brain 
activity in normal arousal states and epileptic seizures. Phys Rev E 65:041924. 

10. Sulfczynski P, Kalitzin S, Lopes Da Silva FH (2004) Dynamics of non-convulsive 
epileptic phenomena modeled by a bistable neuronal network. Neurosci 
126:467-484. 

11. Breakspear M, Roberts JA, Terry JR, Rodrigues S, Mahant N, et al. (2006) A 
unifying explanation of primary generalized seizures through nonlinear brain 
modeling and bifurcation analysis. Cereb Cortex 16:1296—1313. 

12. Lopes da Silva F, Blanes W, Kahtzin SN, ParraJ, Suffczynski P, et al. (2003) 
Epilepsies as Dynamical Diseases of Brain Systems: Basic Models of the 
Transition Between Normal and Epileptic Activity. Epilepsia 44:72-83. 

13. Marten F, Rodrigues S, Benjamin O, Richardson MP, Terry JR (2009) Onset of 
polyspike complexes in a mean-field model of human electroencephalography 
and its application to absence epilepsy. Philos Trans A Math Phys Eng Sci 
28:1145-1161. 

14. Stoeeo A, Lebiere C, .AnelersonJR (2010) Conelitional routing of information to 
the cortex: A model of the basal ganglias role in cognitive coordination. Psychol 
Rev 117:541-574. 

15. Packard MG, Knowlton BJ (2002) Learning and memory functions of the basal 
ganglia. Annu Rev Neurosci 25:563-593. 

16. Groenewegen HJ (2003) The basal ganglia and motor control. Neural Plast 
10:107-120. 

17. Gatev P, Wiehmann T (2009) Interactions between cortical rhythms and spiking 
activity of single basal ganglia neurons in the normal and parkinsonian state. 
Cereb Cortex 19:1330-1344, 



different dynamical states are observed: the saturation state (I), the 
SWD oscillation state (II), the simple oscillation state (III) and the 
low firing state (IV). In each subfigure, the region between two red 
dashed lines denotes the main bidirectional suppression region of 
SWDs, where the double arrow represents that both increasing 
and decreasing Vpii; can inhibit the generation of SWDs. The 
results given in (A)-(H) indicate that the modified model shows 
bistabUity (the simple oscillation state or the saturation state) in 
the large K and strong Vp,; region, and the final dynamics of 
the modified model significantly depend on the initial conditions. 
Note that in aU siiuulations, we set ?o = 80ms, !'(,, = 3.2 mVs, 
Vse = 3.4 mV s, Vre = 1 .6 HiV s, Vjr = — 1 -76 mV s and = 8.0 mV s. 
(TIF) 

Text SI Supporting information code. 

(TXT) 

Text S2 Supporting information code. 

(TXT) 

Acknowledgments 

We sincerely thank Dr. Dan Wu and Dr. Yangsong Zhang for valuable 
comments on early versions of this manuscript. 

Author Contributions 

Conceived and designed the experiments: MC DG TW \'X PX R-WS DY. 
Performed the experiments: MC DG CL. Analyzed the data: MG DG TW 
\\(J YX DY. Contributed reagcnts/malcrials/analysis tools: MC DG. 
Wrote the paper: MG DG YX CL PAVS DY. 



18. Kumar A, Gardanobile S, Rotter S, .Aertsen A (201 1) The role of inhibition in 
generating anel controlling Parkinsons disease related oscillations in the basal 
ganglia. Front Syst Neurosci 5:86. 

19. Paz JT, Deniau JM, Charpier S (2005) Rhythmic bursting in the cortico- 
subthalamo-pallidal network during spontaneous genetically determined spike 
and wave discharges. J Neurosci 25:2092-2101. 

20. Luo C, li Q, Xia Y, Lei X, Xue K, et al. (2012) Resting state basal ganglia 
network in idiopathic generalized epilepsy. Hum Brain Mapp 33:1279-1294. 

21. Deransart C, Vercueil L, Marescaux C, Depaulis A (1998) The role of basal 
ganglia in the control of generalized absence seizures. Epilepsy Res 32:213—223. 

22. Deransart C, Depaulis A (2002) llie control of seizures by the basal ganglia? 
A review of experimental elata. Epileptic Disord. Suppl 3:S61— 72. 

23. Paz JT, Chavez SaiUet S, Deniau JM, Charpier S (2007) .Activity of ventral 
medial thalamic neurons during absence seizures and modulation of cortical 
paroxysms by the nigrothalamic pathway. J Neurosci 27:929-941. 

24. Biraben A, Semah F, Ribeiro MJ, Douaud G, Remy P, et al. (2004) PET 
evidence for a role of the basal ganglia in patients with ring chromosome 20 
epilepsy. Neurology 63:73-77. 

25. Postimia RB, Dagher A (2006) Basal ganglia functional connectivity based on a 
meta-analysis of 126 positron emission tomography anel functional magnetic 
resonance imaging publications. Cereb Cortex 16:1508—1521. 

26. Kase D, Inoue T, Imoto K (2012) Roles of the subthalamic nucleus and 
subthalamic HCN channels in absence seizures. J Neurophysiol 107:393^06. 

27. Haber SN, Calzavara R (2009) The cortico-basat ganglia integrative network: 
the role of the thalamus. Brain Res Bull 78:69-74. 

28. Gulcebi MI, Ketenci S, Linke R, Hacioglu H, Yanali H, et al. (2012) 
Topographical connections of the substantia nigra pars reticttlata to higher-order 
thalamic nuclei in the rat. Brain Res Bull 87:312—318. 

29. Di Chiara G, Porceddu ML, MoreUi M, Mulas ML, Gessa GL (1979) Evidence 
for a (iABAergie projection from the substantia nigra to the ventromedial 
thalamus anel to the superior coUiculus of the rat. Brain Res 176:273 -284. 

30. van .\lbada SJ, Robinson PA (2009) Mean-field modeling of the basal ganglia- 
thalamoeortieal svstem. I: Firing rates in healthy and parkinsonian states. J Theor 
Biol 257:642-663. 

31. van Albada SJ, Gray RT, Drysdale PM, Robinson PA (2009) Mean-field 
modeling of the basal ganglia-thalamocortical system. II: Dynamics of 
parkinsonian oscillations. J Theor Biol 257:664—688. 

32. Bar-Gad I, Morris G, Bergman H. (2003) Information processing, dimension- 
ality reduction and reinforcement learning in the basal ganglia. Prog Neurobiol 
71:439-473. 

33. Robinson PA, Rennie CJ, Wright JJ (1997) Propagation and stability of waves of 
electrical activity in the cerebral cortex. Phys Rev E 56:826-840. 



PLOS Computational Biology | www.ploscompbiol.org 



16 



March 2014 | Volume 10 | Issue 3 | el 003495 



Control of Absence Seizures by the Basal Ganglia 



34. Robinson PA, Rcnnic GJ, Wright JJ, Bahramali H, Gordon E, ct al. (2001) 
Prediction of electroenccphalographic spectra from neurophysiology. Phys Rev E 
63:021903. 

35. Robinson PA, Rennie CJ, Wright JJ, Bourke PD (1998) Steady states and global 
dynamics of electrical activity in the cerebral cortex. Phys Rev E 58:3557—3571. 

36. Robinson PA, Rennie CJ, Rowe DL, O'Connor SC (2003) Estimation of 
multiscale neurophysiologic parameters by electroenccphalographic means. 
Hum Brain Mapp 23:53-72. 

37. Dcstcxhc A. Srjnowski T] (2001) Thalamocortical assemblies. Oxford: Oxford 
Unix. Press. 

38. Jirsa VK, Hakrn H (1996) field theory of electromagnetic brain activity. Phys 
Rev Lett 77:960-963. 

39. Goodfellow M, Schindler K, Baier G (2011) Intermittent spikc-wavc dynamics 
in a heterogeneous, spatially extended neural mass model. Neuroimage 20 1 1 
55:920-932. 

40. Butcher JC (2003) Numerical Methods for Ordinary Differential Equations, 
New York: John Wiley & Sons. 7-1 1 p. 

41. Ermcntrout (2002) Simulating, Analyzing, and Animating Dynamical Systems. 
PhUadelphia PA: SIAM. 

42. Hosfbrd DA, Clark S, Cao Z, Wilson WAJr, lin FH (1992) The role of GABAb 
receptor activation in absence seizures of lethargic (Ih/lh) mice. Science 
257:398-401. 

43. Dcstcxhc A (2008) Gortico-thalamic feedback: a key to explain absence seizures. 
In: Soltesz I, Staley K, eds. Computational Neuroscicncc in Epilepsy. Elsevier. 
184-214 p. 

44. Raastad M, Storm JF, Andersen P (1992) Putative Single Quantum and Single 
Fibre Excitatory Postsynaptic Currents Show Similar Amplitude Range and 
Variability in Rat Hippocampal Slices. EurJ Neurosci 4:113-117. 

45. Roscnmund G, Clements JD, Wcstbrook GL (1993) Nonuniform probability of 
glutamate release at a hippocampal synapse. Science 262:754-757. 

46. Allen C, Stevens GF (1994) An evaluation of causes for unreliability of synaptic 
transmission. Proc Natl Acad Sci USA 91:10380-10383. 

47. Goldman MS, Maldonado P, Abbott LF (2002) Redundancy rrduedon and 
sustained firing with stochastic depressing synapses. J Neurosci 22:584 -591. 

48. Guo DQj Li GG (2011) Signal propagation in feedforward neuronal networks 
with unreliable synapses, J Comput Neurosci 30:567—587. 

49. Guo DQ, Li CG (2012) Stochastic resonance in Hodgkin-Huxley neuron 
induced by unreliable synaptic transmission, J Theor Biol 308:105-114. 

50. Perry MS, Bailey LJ, Kotecha AC, Malik SI, Hernandez AW (2012) 
Amantadine for the treatment of refractory absence seizures in children. Pediatr 
Neurol 46:243-245. 



51. Bouma PA, Westcndorp RG, van Dijk JG, Peters AG, Brouwer OF (1996) 
The outcome of absence epilepsy: a meta-analysis. Neurology 47:802—808. 

52. Vereucil L, Bcnazzouz A, Deransart C, Bressand K, Marescaux C, et al. (1998) 
High-frequency stimulation of the subthalamic nucleus suppresses absence 
seizures in the rat: comparison with neurotoxic lesions. Epilepsy Res 31:39^6. 

53. Krampfl K, Maljevic S, Cossette P, Ziegler E, Rouleau GA, et al. (2005) 
Molecular analysis of the A322D mutation in the GABA receptor alpha-subunit 
causing juvenile myoclonic epilepsy, EurJ Neurosci 22:10 20. 

54. Hie A, Raimondo JV, Akerman GJ (2012) Adenosine release during seizures 
attenuates (iABAA receptor-mediated depolarization. J Neurosci 32:5321-5332. 

55. Frdhlich F, Scjnowski IJ, Bazhenov M (2010) Network bistability mediates 
spontaneous transitions between normal and pathological brain states. J Neurosci 
30:10734-10743. 

56. Freyer F, Roberts JA, Becker R, Robinson PA, Ritter P, et al. (201 1) Biophysical 
mechanisms of multistability in resting-state cortical rhythms. J Neurosci 

31:6353-6361. 

57. Volman V, Pcrc M, Bazhenov M (2011) (jap junctions and epileptic seizures- 
Two sides of the same coin?, PLoS ONE 6:e20572. 

58. Netoff TI, Glewley R, Arno S, Keck T, White JA (2004) Epilepsy in small-world 
networks. J Neurosci 24:8075-8083. 

59. Ponten SC, Bartolomei F, Stam CJ (2007) Small-world networks Euid epilepsy: 
graph theoretical analysis of intracerebrally recorded mesial temporal lobe 
seizures. Clin Neurophysiol 118:918-927. 

60. Nail-Boucherie K, Le-Pham BT, Marescaux C, Depaulis A (2002) Suppression 
of absence seizures by electrical and pharmacological activation of the caudal 
superior colliculus in a genetic model of absence epilepsy in the rat. Exp Neurol 
177:503-514. 

61. Wang Chen G, Perc M (2011) Synchronous bursts on scale-free neuronal 
networks with attractive and repulsive coupling. PLoS One 6:el5851. 

62. Wilson HR, Cowan JD (1972) Excitatory and inhibitory interactions in localized 
populations of model neurons. BiophysJ 12: 1—24. 

63. Wilson HR, Cowan JD (1973) A mathematical theory of the functional dynamics 
of cortical and thalamic nervous tissue. Kybernctik 13: 55-80. 

64. Taylor PN, GoodfeUow M, Wang Y, Baier G (2013) Towards a large-scale 
model of patient-specific epileptic spike-wave dischju-ges, Biol Cybem 107:83- 
94. 

65. Friston KJ, Harrison L, Penny W (2003) Dynamic causal modelling. Neuro- 
image 19:1273-1302. 

66. Marreiros AC, Kiebel SJ, Friston KJ (2010) A dynamic causal model study of 
neuronal population dyn£unics. Neuroimage 51:91-101. 

67. Pinotsis DA, Moran RJ, Friston KJ (2012) Dynamic causal modeling with neural 
fields. Neuroimage 59:1261-1274. 



PLOS Computational Biology | www.ploscompbiol.org 



17 



March 2014 | Volume 10 | Issue 3 | e1003495 



